Systems with Multiplicative Noise: Critical Behavior from KPZ Equation and 

Numerics 
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We show that certain critical exponents of systems with multiplicative noise can be obtained from 
exponents of the KPZ equation. Numerical simulations in Id confirm this prediction, and yield other 
exponents of the multiplicative noise problem. The numerics also verify an earlier prediction of the 
divergence of the susceptibility over an entire range of control parameter values, and show that the 
exponent governing the divergence in this range varies continuously with control parameter. 
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Langevin equations - first order partial differential 
equations in time containing Gaussian random noise 
terms - capture the macroscopic physics of many clas- 
sical, stochastic, many-body systems JjJ. In the most 
common situation, one that includes much of equilib- 
rium statistical mechanics, the noise amplitude is sim- 
ply a constant. There are however, important classes of 
problems in which the noise amplitude is proportional to 
a positive power, a, of the field variable itself. The well- 
known case a = 1/2 describes the physics of directed 
percolation and its many and diverse realizations Q in, 
e.g., models of catalysis, the spread of epidemics or forest 
fires, and propagating star formation. Other systems in 
which the dominant source of noise is external, among 
them certain chemical reactions and problems in quan- 
tum optics, are characterized by a = 1 ||. These so- 
called "multiplicative-noise" systems (MNS) || are the 
subject of this paper. 

The generic problem with a single-component field in 
this class is described by the Langevin equation 



dn{x,t) / dt = V 2 n — rn 



+ nrj 



(1) 



where r and u are parameters, and 77 is a Gaussian 
noise variable with correlations < r)(x,t)r)(x' ,<') >= 
DS(x — x')8{t — t'), for some noise strength D. We con- 
sider both quadratic and cubic nonlinearities, p = and 
p = 1, respectively. 

In a previous paper j5j, we analyzed the phase struc- 
ture of this model, which has two phases, an "active" 
phase with < n > > 0, and an "absorbing" phase with 
n{x) = for all x, that occurs for sufficiently large r. 
These phases are separated by a critical point at r = r c , 
where r c — in mean-field theory. (Note that the van- 
ishing of the noise amplitude with n makes the right side 
of Eq. (Q) vanish when n(x) vanishes, thereby causing 
all dynamics to cease and making the absorbing configu- 
ration n(x) = a potentially stable phase.) We showed 
that model (|l|) had a critical dimension d c — 2, below 
which the transition is governed by an analytically in- 
accessible strong-coupling fixed point. For d > 2, on 
the other hand, the transition is governed by the "weak- 
coupling" Gaussian fixed point with mean-field expo- 



nents for D less than a critical value D c . For D > D c , 
however, the strong-coupling fixed point is the stable one. 
A multicritical point occurs at D = D c . 

In this paper we argue that the critical behavior of 
MNS is actually governed by the fixed point of the 
Kardar-Parisi-Zhang (KPZ) model of growing interfaces 
||. This mapping is consistent with the phase diagram 
proposed in ref. H, and with known exponents of the 
weak-coupling and multicritical fixed points. It also al- 
lows us to express the dynamical exponent z and the cor- 
relation length exponent v of the strong-coupling transi- 
tion in terms of the KPZ exponents. These exponents 
are found to be independent of the degree, 2 + p, of 
the nonlinearity. A lower bound of unity for the order 
parameter exponent (3 in the case p = also emerges. 
We confirm these predictions by calculating numerically 
in Id the four independent exponents characterizing the 
strong-coupling transition. We also confirm the rather 
striking prediction in ref. || of an entire domain of r 
values, encompassing the critical value r c , in which the 
susceptibility of the system diverges. The numerics show 
that, like in the exactly solvable single variable (Od) prob- 
lem, this region of infinite susceptibility extends to both 
sides of the critical point, and is controlled by a fixed line 
with continuously varying critical exponents Q . 

To establish the connection between model (1) and the 
KPZ theory, note that the field n(x, t) in (1) will remain 
positive if n(x, 0) > for all x. In this case one can per- 
form the Hopf-Cole Rj change of variable n(x, t) = e h ^ x,t ^ , 
producing the equation @j 



dh/dt = -r + W 2 h+(Wh) 2 



V 



(2) 



Aside from the u term this is precisely the KPZ equa- 
tion (wherein the standard KPZ nonlinearity (V/1) 2 has 
coefficient unity). Note however that either in the ab- 
sorbing phase or at the critical point, the steady-state 
value of n is zero, whereupon the steady-state value of h 
is —00. Thus the u term vanishes in steady state, leaving 
one with precisely the KPZ theory. 

Recall H that the phase diagram of the KPZ equation 
consists of a unique strong coupling phase for d < d c = 2, 
and both weak and strong coupling phases separated by 
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a multicritical point for d > 2. The weak and strong- 
coupling regimes occur for noise strengths D that are re- 
spectively smaller and larger than a critical value. This 
phase diagram is thus encouragingly similar to that of 
the MNS. Also, the dynamical critical exponent z is 2 for 
both MNS H and KPZ g along the line of weak-coupling 
transitions and at the multicritical point, for d > 2. 

Let us now consider the critical exponents for the 
strong-coupling transition. The dynamical critical ex- 
ponent z for the MNS can be computed from the steady- 
state behavior of response functions right at the critical 
point, and hence is identical to the value of z in the KPZ 
theory. In particular then, z = 3/2 for d = 1 ||. 

To argue that other exponents of the MNS can also 
be obtained from the KPZ equation requires us to con- 
sider the active state of model (|l|). For r slightly less 
than r c , < < n > < 1 in steady state, imply- 
ing that ho =< h > is very large and negative. Writing 
h(x,t) =< h > +Sh(x, t), one obtains an equation for Sh 
that is identical to the KPZ equation except for the extra 
nonlinear term — u'e^ 1+f ^ sh , where v! = ue^ 1+p ^ h °. Not- 
ing that the leading nontrivial term in the expansion of 
this nonlinearity in powers of Sh is the linear "mass" term 
— u'(l + p)Sh, we infer that the main effect of — u'e^ 1+p ^ Sh 
is to produce a finite correlation length £ at which the 
power law correlations of the KPZ equation are cut off 
and replaced by exponential behavior. One concludes 
that £ must be the correlation length of the correspond- 
ing MNS; its divergence as r — > r c is governed by the 
critical exponent v. In the critical region, i.e., on length 
scales |x| << £, the u' term is negligible, so critical cor- 
relations can be computed from the KPZ equation. 

To calculate v from KPZ, take the expectation value 
of the Sh equation, recall that < Sh >= 0, and write the 
extra nonlinear term as —u < n 1+p >, to obtain 

- r+ < (Vh) 2 > -u < n 1+p >= . (3) 

At the critical point, r = r c and n = 0, so — r c + < 
(V/i) 2 > c — . Subtracting this equation from (||) yields 
SW = —5r + u < n 1+p >, where 5r = r c ~r, and SW =< 
(Vh) 2 > - < (Vh) 2 > c . ^From the standard scaling of 

the KPZ equation [§, one has SW Cf 2 ^- 1 ), where 

C is a positive constant, and x the roughness exponent 
of a KPZ interface; i.e., < {h(x, t) - h(0, t)) 2 >~ x 2 * for 
KPZ in steady-state. Since C > 0, Sr > in the active 
phase, and < n 1+p >~ (Sr) 01 +i> > 0, the equation for £ 
has a solution for small Sr only if the exponent (3\+ p is 
greater than unity. For the quadratic nonlinearity, p = 0, 
so < n l+p > is the order parameter < n >, whereupon 
this constraint places the nontrivial bound j3 > 1 on the 
order parameter exponent = f3\. For any value of p, 
one then obtains the result £ ~ Sr~ v with v = 1/(2 — 2^)- 
(Note that in cases where the KPZ interface is smooth, 
X = in this formula.) 

We now describe the numerical simulation of ([l]) in Id. 
We discretize the continuum equation as: 

rii(t + At) = m(t) + At[-rm(t) - uni(tf +p 



+ (l/Ax 2 ){n i+1 {t)+n i - 1 (t)-2n i {t))] 

+ VDAtm(t)rh(t) (i = l,2,-..,JV), (4) 

where Ax is the lattice spacing, At is the time step and 
r]i(t) is a Gaussian random number with unit standard 
deviation. By choosing the argument of the field mul- 
tiplying the noise as in Eq. (0), we are using the Ito 
interpretation [Q of the stochastic differential equation. 
We set Ax = 1 , At = 0.02, and use periodic boundary 
conditions with the system size L = NAx. We also fix 
the noise amplitude \f~D — 4, and vary the linear coeffi- 
cient r as the only control parameter. In the following, 
we present our numerical results for p = 1 fl(][| . 

In numerical simulation, due to the finite time step 
At, the property of (1) that n(x,t) > if n(x, 0) > 
for all x is lost. However, it is easy to fix this prob- 
lem by setting rii(t) to zero if its value becomes neg- 
ative under Eq. (El). The effect of this modification 
is easily estimated. For a single time step, we can ne- 
glect the second term on the RHS of (^) because it is 
of higher order in At than the noise term. Then set- 
ting negative values of rii (t) to is equivalent to truncat- 
ing the probability distribution of rji{t) so that its mini- 
mum is set by rj m i n = — l/V-DAi, and replacing all the 
V < Vmin by r\ m in- This means, however, that the mean 
of rji(t) is no longer zero, and thus a deterministic term 
proportional to Uj(t) is generated. The resulting effec- 
tive linear coefficient r e ff can be roughly estimated as: 
r eff = r-^/D/2^At J*™ n (r) min -r)) exp(-?7 2 /2)^. The 
effective strength of the noise is also changed because of 
the truncation. These changes in parameters should not, 
however, alter the universality class of the transition [ pd| . 

The first step in studying the critical behavior of 
MNS numerically is to locate the critical point. Starting 
with uniform initial conditions, and letting the equation 
evolve long enough to reach steady state, we compute 
M =< rii(t) > for different values of r, where < > 
denotes both spatial and temporal averages, as well as 
averages over different independent runs. We average 
over between 2* 10 5 and 6* 10 6 time steps, and up to 100 
independent runs, depending on the system size. To han- 
dle finite-size effects, we studied different system sizes: 
N = 100,200,400,1000. In fig. 1(a), we show the de- 
pendence on the inverse system size 1 /N of the critical 
point, r c (N), defined by M first becoming numerically 
indistinguishable from zero in every run. Extrapolating 
the fitted line [Q to N — oo determines a critical value 
r c w —2.18. Fig. 1(b) shows M vs. Sr = r c — r on a 
log-log plot for N=400. The best fit to M - (Sr) yields 
||| /3= 1.70 ±0.05. 

Also depicted in fig. 1(b) are the higher order moments 
of the n field: M m —< n m >, with m = 2, 3. Evidently, 
M2 and M3 also have a power law dependence on Sr: 
Mm ~ (8r)P m with /3 2 and fa being equal to within 
our numerical accuracy. It is not surprising that there 
is anomalous scaling, i.e., 02 ^ 2/3, (3 3 ^ 3/3, and so on, 
since the strong coupling fixed point is non-Gaussian, but 
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it remains to be seen whether (3 m is indeed independent of 
to. It is interesting to note that for the zero-dimensional 
(single-variable) case, where an exact solution is avail- 
able BEl all the moments do scale with exactly the 
same exponent: < n m >= 2 m / 2 T [ m +^/ D ] /T{-r/2D) , 
i.e., (3 m = 1 for all to. 

To determine the other critical exponents, we need 
to calculate the two-point function C{x,t) =< nix + 
xo,t + to)n(xo,to) > in the active phase. We 
have computed both C(x, 0) and C(0, f), for — r = 
2.8, 2.6, 2.5, 2.4, 2.34, 2.3, 2.2, 2.215. We used N=512, and 
averaged over T — 10 6 At time steps in steady state, and 
over up to 100 independent runs. 

We summarize the results in fig. 2. In fig. 2(a,b), 
C(x,0) and C(0, i) are plotted for different values of 
r. It is clear from these figures that the scaling regime 
becomes bigger as we approach the critical point, and 
that the amplitudes of the correlation functions become 
smaller, consistent with C(x,t) vanishing, as it must, in 
the absorbing phase. In fig. 2(c), we plot the abso- 
lute value of the connected space and time correlation 
function, C c (x,t) — C(x,t) — M 2 , at the largest value 
of r, viz., -2.215. The power law decay of the correla- 
tion functions in the scaling regime can be characterized 
by C c (0,i) = A t (r)t~ at and C c (x,0) = A x {r)x~ a * with 
the exponents at = 1.08 ± 0.04 and a x = 1.65 ± 0.07. 
In conventional notation, this corresponds to exponents 
7] = 2-d+a x = 2.65±0.07, and z = a x /a t = 1.53±0.07. 
In fig. 2(d), we plot the dependence of the amplitudes 
of the spatial and temporal correlation functions, A x (r) 
and A t (r), on Sr in the scaling regime; these can be fitted 
by A x . t (r) ~ Sr A with the exponent A = 1.7 ± 0.07. 

^From these measured values of the four indepen- 
dent exponents (3, rj, z, and A, we obtain the corre- 
lation length exponent v from the scaling relation || 
v= (2(3- A)/(d-2 + r)) = 1.03 ±0.05. These values of v 
and z are in excellent agreement with the values, 1 and 
3/2 respectively, following from our argument that the 
critical behavior of MNS is controlled by the KPZ equa- 
tion; /?2 is also > 1, consistent with the bound derived 
from KPZ. To further check the accuracy of various scal- 
ing exponents, we have also measured the decay of the 
average density right at the critical point, starting from 
a homogeneous initial condition: M ~ t~ e . Using scal- 
ing arguments, it is easy to express in terms of other 
exponents: 8 — (3/(vz). From the numerical values of (3, 
v, and z, we predict 8 to be 1.079. In fig. 3, we plot 
M versus time at r = r c , and the exponent 9 thus mea- 
sured is 9 = 1.1 ± 0.05, in excellent agreement with the 
scaling prediction. In the same figure, we have included 
the behavior of the higher order moments M2$,i(t) at 
the critical point. These plots strongly suggest that the 
exponents 8 m = (3 m j(yz) for these higher moments are 
all equal to 8, in agreement with the static measurement. 

Another important characterization of any phase tran- 
sition is the response function. For equilibrium sys- 
tems, the fluctuation-dissipation theorem ensures that 



the response function is directly related to the correlation 
function of the system, so the susceptibility probes the 
system's internal fluctuations and correlations directly. 
However, the fluctuation-dissipation theorem does not 
generally hold for nonequilibrium systems, where the re- 
sponse and correlation functions can differ significantly. 
Indeed, one of the striking features of MNS is that the 
susceptibility is predicted |5[ to diverge over a finite range 
of control parameters, while the correlation function has, 
more conventionally, a unique critical point. 

To test this prediction numerically in Id, we add a con- 
stant source term <j> on the RHS of Eq. (||) , and calculate 
the average density M(</>, r) for small values of </> and for 
each value of r. In figure 4(a), we show the dependence of 
M((j),r) on <p for different values of r around the critical 
point r c . We can fit the curves for small 4> by M (</>, r) = 
M (0, r) + B(r)<jP^ for small <f>, where M(0, r) is the or- 
der parameter in the absence of any external source, and 
B(r) is a function of r. The static susceptibility is defined 
as x = lim^o dM/dip — lim^o -B(r , )7(r)0 7 ( r )~ 1 , so the 
susceptibility diverges when j(r) < 1. According to fig. 
4(b), where we show the dependence on r of the sus- 
ceptibility exponent, 7(7") — 1, x diverges not just in the 
absorbing phase, but in the whole region 3.2 > r > —5.9 
surrounding the critical point r c = —2.18. Both the con- 
tinuous variation of this exponent with r, and the diver- 
gence of x m both the active and absorbing phases near 
r c also occur in the exactly solvable Od problem [H. 

Another quantitative prediction of our previous paper 
p] is that the region of diverging x should terminate in 
the absorbing phase when the bare mass vanishes, i.e., 
at r = in the Ito representation. However, because of 
the shift in r due to our numerical algorithm, the lower 
boundary of this region is shifted to r > 0. A clearer 
way to test the prediction is therefore to make the Hopf- 
Cole transformation n = e h first, and then simulate the 
equation for h, noting that the external field will have 
the form 4>e~ h . Because the Hopf-Cole relation automat- 
ically guarantees that the n field is positive definite, there 
is no need to alter the numerical algorithm, so 7(r = 0) 
should equal 1. We have verified this in a simulation. 

We thank A. Pikovsky, P. Grassberger, T. Hwa, S. 
Chen, and N.-Z. Cao for helpful discussions, and T. Hwa 
for bringing refs. || to our attention. 
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FIG. 1. (a)r c (N) vs. inverse system size 100/N; (b)< n m > 
with m = 1,2,3 vs. control parameter 8r = r c — r; dashed 
line with slope 1.7 gives the best fit to the data. 

FIG. 2. (a) and (b) show the spatial and temporal corre- 
lation functions with different values of r < r c ; (c) connected 
correlation functions at r = —2.215; dashed lines are fits with 
slopes 1.08 and 1.65; (d) amplitudes of the correlation func- 
tions vs. 5r; dashed line has slope 1.7. 

FIG. 3. Power-law decay in time of various moments of the 
order parameter at the critical point r = r c . 

FIG. 4. (a) Order parameter vs. external field at values of 
r near r c ; dashed lines are fits to M(<f>,r) — A/(0, r) ~ 
as explained in text; (b) susceptibility exponent vs. r; sus- 
ceptibility diverges for a range of r values. 
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